## Replication code for Section 3 of the Supplmentary Information for:
##  Abrajano, Elmendorf, and Quinn.
##  "Labels vs. Pictures: Treatment Mode Effects in Experiments About
##   Discrimination." Political Analysis. 
##
## Marisa Abrajano, Christopher Elmendorf, and Kevin Quinn
## 9/24/2017
##


library(lmtest)
library(sandwich)
library(msm)
library(ggplot2)
library(car)


## function that does clustered SEs
## (from Hainmueller, Hopkins, & Yamomoto)
vcovCluster <- function(model, cluster){
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}


## load the long (stacked) data
load("./PhotosLong-ALL.Rda")
load("./WordsLong-ALL.Rda")

## convert choice variable into numeric
mydata.w.long$choice <- as.numeric(mydata.w.long$choice == "CandA")
mydata.p.long$choice <- as.numeric(mydata.p.long$choice == "CandA")

## convert matchup.time to factor
mydata.w.long$matchup.time <- factor(mydata.w.long$matchup.time,
                                     levels=c("FIRST", "SECOND", "THIRD",
                                              "FOURTH", "FIFTH", "SIXTH"))
mydata.p.long$matchup.time <- factor(mydata.p.long$matchup.time,
                                     levels=c("FIRST", "SECOND", "THIRD",
                                              "FOURTH", "FIFTH", "SIXTH"))


###########################################################################
## Hypothesis 1 (and 4 and 5)

## estimate tau (effect of text ethnicity on choice for full sample)
mydata.sub <- na.omit(mydata.w.long)
lm.tau <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.tau <- coeftest(lm.tau, vcov=vcovCluster(lm.tau,
                                             cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.tau))
wald.tau <- linearHypothesis(lm.tau, vcov=vcovCluster(lm.tau,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])
                              



## estimate nu (effect of photo ethnicity on choice for full sample)
mydata.sub <- na.omit(mydata.p.long)
lm.nu <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.nu <- coeftest(lm.nu, vcov=vcovCluster(lm.nu,
                                           cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.nu))
wald.nu <- linearHypothesis(lm.nu, vcov=vcovCluster(lm.nu,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])





###########################################################################
## Hypothesis 2

## estimate tau_I (effect of text ethnicity on choice for high int. mot. )
mydata.sub <- mydata.w.long[mydata.w.long$int.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.tauI <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.tauI <- coeftest(lm.tauI, vcov=vcovCluster(lm.tauI,
                                             cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.tauI))
wald.tauI <- linearHypothesis(lm.tauI, vcov=vcovCluster(lm.tauI,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])
V <- vcovCluster(lm.tauI, cluster=mydata.sub$respondent.ID)
beta.hat <- coef(lm.tauI)
save(V, file="./V-tauI.Rda")
save(beta.hat, file="./beta-tauI.Rda")



## estimate tau_notI (effect of text ethnicity on choice for low int. mot. )
mydata.sub <- mydata.w.long[mydata.w.long$int.mot.high==FALSE,]
mydata.sub <- na.omit(mydata.sub)
lm.taunotI <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.taunotI <- coeftest(lm.taunotI, vcov=vcovCluster(lm.taunotI,
                                             cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.taunotI))
wald.taunotI <- linearHypothesis(lm.taunotI, vcov=vcovCluster(lm.taunotI,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])


## estimate nu_I (effect of photo ethnicity on choice for high int. mot. )
mydata.sub <- mydata.p.long[mydata.p.long$int.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.nuI <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.nuI <- coeftest(lm.nuI, vcov=vcovCluster(lm.nuI,
                                             cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.nuI))
wald.nuI <- linearHypothesis(lm.nuI, vcov=vcovCluster(lm.nuI,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])

V <- vcovCluster(lm.nuI, cluster=mydata.sub$respondent.ID)
beta.hat <- coef(lm.nuI)
save(V, file="./V-nuI.Rda")
save(beta.hat, file="./beta-nuI.Rda")


## estimate nu_notI (effect of photo ethnicity on choice for low int. mot. )
mydata.sub <- mydata.p.long[mydata.p.long$int.mot.high==FALSE,]
mydata.sub <- na.omit(mydata.sub)
lm.nunotI <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.nunotI <- coeftest(lm.nunotI, vcov=vcovCluster(lm.nunotI,
                                             cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.nunotI))
wald.nunotI <- linearHypothesis(lm.nunotI, vcov=vcovCluster(lm.nunotI,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])





###########################################################################
## Hypothesis 3

## estimate tau_E (effect of text ethnicity on choice for high ext. mot. )
mydata.sub <- mydata.w.long[mydata.w.long$int.mot.high==FALSE & mydata.w.long$ext.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.tauE <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.tauE <- coeftest(lm.tauE, vcov=vcovCluster(lm.tauE,
                                             cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.tauE))
wald.tauE <- linearHypothesis(lm.tauE, vcov=vcovCluster(lm.tauE,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])


## estimate nu_E (effect of photo ethnicity on choice for high ext. mot. )
mydata.sub <- mydata.p.long[mydata.p.long$int.mot.high==FALSE & mydata.p.long$ext.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.nuE <- lm(choice ~ candA.ethnicity*matchup.time, data=mydata.sub)
tab.nuE <- coeftest(lm.nuE, vcov=vcovCluster(lm.nuE,
                                             cluster=mydata.sub$respondent.ID))
beta.hat <- names(coef(lm.nuE))
wald.nuE <- linearHypothesis(lm.nuE, vcov=vcovCluster(lm.nuE,
                                                       cluster=mydata.sub$respondent.ID),
                              hypothesis.matrix=beta.hat[grep(":", beta.hat)])







load("V-nuI.Rda")
load("beta-nuI.Rda")

rounds <- c("First", "Second", "Third", "Fourth", "Fifth", "Sixth")
xvals <- 1:6

nu.hat <- rep(NA, 6)
nu.hat[1] <- beta.hat[2]
nu.hat[2] <- beta.hat[2] + beta.hat[8]
nu.hat[3] <- beta.hat[2] + beta.hat[9]
nu.hat[4] <- beta.hat[2] + beta.hat[10]
nu.hat[5] <- beta.hat[2] + beta.hat[11]
nu.hat[6] <- beta.hat[2] + beta.hat[12]

nu.se <- rep(NA, 6)
nu.se[1] <- V[2,2]
nu.se[2] <- V[2,2] + V[8,8] + 2*V[2,8]
nu.se[3] <- V[2,2] + V[9,9] + 2*V[2,9]
nu.se[4] <- V[2,2] + V[10,10] + 2*V[2,10]
nu.se[5] <- V[2,2] + V[11,11] + 2*V[2,11]
nu.se[6] <- V[2,2] + V[12,12] + 2*V[2,12]

nu.se <- sqrt(nu.se)

pdf("nuI-tauI-bytime.pdf", height=5.5, width=12)
par(mfrow=c(1,2), mar=c(5, 6, 2, 4)+.1, lend=2)
plot(xvals, xvals, ylim=c(-.2, 0.4), type="n",
     xlab="Round", ylab="Latino Candidate Effect Size",
     xaxt="n", main=expression(nu[I]))
abline(h=seq(from=-1, to=1, by=0.05), col="gray96")
abline(h=seq(from=-1, to=1, by=0.1), col="gray94")
abline(h=seq(from=-1, to=1, by=0.2), col="gray93")
abline(h=0, col="gray88")

axis(side=1, at=xvals, labels=rounds, cex.axis=.8)

nu.col <- rgb(1, 0.4098, 0.1784)
points(xvals, nu.hat, pch=16, col=nu.col)
segments(x0=xvals, x1=xvals, y0=nu.hat-1.96*nu.se,
         y1=nu.hat+1.96*nu.se, col=nu.col, lwd=2.5, lend=2)





load("V-tauI.Rda")
load("beta-tauI.Rda")

rounds <- c("First", "Second", "Third", "Fourth", "Fifth", "Sixth")
xvals <- 1:6

tau.hat <- rep(NA, 6)
tau.hat[1] <- beta.hat[2]
tau.hat[2] <- beta.hat[2] + beta.hat[8]
tau.hat[3] <- beta.hat[2] + beta.hat[9]
tau.hat[4] <- beta.hat[2] + beta.hat[10]
tau.hat[5] <- beta.hat[2] + beta.hat[11]
tau.hat[6] <- beta.hat[2] + beta.hat[12]

tau.se <- rep(NA, 6)
tau.se[1] <- V[2,2]
tau.se[2] <- V[2,2] + V[8,8] + 2*V[2,8]
tau.se[3] <- V[2,2] + V[9,9] + 2*V[2,9]
tau.se[4] <- V[2,2] + V[10,10] + 2*V[2,10]
tau.se[5] <- V[2,2] + V[11,11] + 2*V[2,11]
tau.se[6] <- V[2,2] + V[12,12] + 2*V[2,12]

tau.se <- sqrt(tau.se)

plot(xvals, xvals, ylim=c(-.2, 0.4), type="n",
     xlab="Round", ylab="Latino Candidate Effect Size",
     xaxt="n", main=expression(tau[I]))
abline(h=seq(from=-1, to=1, by=0.05), col="gray96")
abline(h=seq(from=-1, to=1, by=0.1), col="gray94")
abline(h=seq(from=-1, to=1, by=0.2), col="gray93")
abline(h=0, col="gray88")

axis(side=1, at=xvals, labels=rounds, cex.axis=.8)

nu.col <- rgb(1, 0.4098, 0.1784)
tau.col <- rgb(0.0627, 0.3059, 0.54511)

points(xvals, tau.hat, pch=16, col=tau.col)
segments(x0=xvals, x1=xvals, y0=tau.hat-1.96*tau.se,
         y1=tau.hat+1.96*tau.se, col=tau.col, lwd=2.5, lend=2)



dev.off()









